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1.    Introduction 

The  central  issue  in  approximating  a  partial  differential  equation  by 
a  finite  difference  equation  is  the  degree  to  which  the  difference  equation 
solution  agrees  with  the  solution  to  the  differential  equation.   This 
agreement  can  be  considered  in  both  its  quantitative  (e.g.,  relative 
error),  or  qualitative  aspects,  e.g.,  behavior  of  transients,  propagation 
of  fronts,  etc. 

In  this  paper  we  investigate  some  of  these  aspects  for  a  prototype 
dispersive  wave  (hyperbolic)  equation.   This  type  of  problem  is 
especially  important  in  numerical  weather  prediction,  since  for  "fast" 
time  scales,  the  meteorological  primitive  equations  are  dispersive.   We 
consider  both  the  differential  and  quasi-discrete  case,  where  second  order 
centered  spatial  differences  are  used.   In  both  cases  solutions  are 
obtained  by  Fourier  transform  methods,  with  steady  state  solutions 
extracted  directly  from  the  transforms  and  inverted  in  closed  form,  and 
the  asymptotic  behavior  of  the  transients  determined  by  the  method  of 
stationary  phase. 

Our  analysis  shows  that  discretization  of  spatial  derivatives  has 
two  major  effects: 

(1)  The  initial  conditions  of  the  differential  problem  contribute 
to  the  steady  state  solution  in  a  manner  that  decays  exponentially  with 
distance.   This  qualitative  effect  is  retained  in  the  quasi-discrete 
formulation,  but  the  rate  of  decay  is  decreased. 

(2)  The  discretization  introduces  additional  transients  beyond 
those  encountered  in  the  differential  case.   These  new  transients  are  at 


least  of  the  same  magnitude  of  decay  as  the  differential  transients,  and, 
in  some  instances,  they  are  more  persistent. 


2.   The  Differential  One-Dimensional  Adjustment  Process 

One  of  the  simplest  models  of  dispersive  waves  is  the  linearized  one- 
dimensional  shallow  water  equations  with  no  mean  flow,  in  an  infinite 
region: 

|»-*  +  g*.o  (i) 

|*  +  fu  -  0  (2) 

where  u  is  the  perturbation  velocity  in  the  x  direction,   v  is  the 
perturbation  velocity  normal  to  the  x  direction,   H  and  h  are  the  mean 
and  perturbed  heights  of  the  free  surface,  respectively,  and  g  >  0  and 
f  >  0  are  gravitational  and  Coriolis  parameters,  respectively.   This  model 
is  especially  important  in  the  study  of  the  meteorological  problem  of 
geostrophic  adjustment,  and  has  been  studied  in  some  detail  by  Rossby  [1], 
Cahn  [2],  Blumen  [3]  and  Winninghoff  [4].   In  their  papers,  the  model  has 
been  studied  by  eliminating  between  the  equations  to  arrive  at: 

2  2 

9  u  ,  c2        9  u   A  ,/N 

— -  +  f   u  -  gH  — *  =  0  ,  (4) 

3t  3x 

then  solving  (4)  by  a  Fourier  Transform  approach.   After  solving  (4), 
solutions  for  h  and  v  are  obtained  by  substitution  into  (2)  and  (3), 
although  closed- form  solutions  are  not  produced  in  some  of  the  papers. 
Note  the  dispersive  character  of  (4)  is  clearly  seen  by  assuming  a  wave 
solution: 


t      .N    .   i(kx-vt) 
u(x,t)  =  A  e  , 

which  leads  immediately  to: 

v2  =   f2  +  k2gH  -   f2(l+A2k2),      X   =   /gll/f    .  (5) 

We  now  derive  an  alternative  means  of  solving  (l)-(3)  which  is 
superior  in  that  it  does  not  require  elimination,  it  produces   u,  v  and  h 
without  back  substitution,  it  yields  interesting  insights  into  the 
transient  and  steady  state  behavior  of  the  solutions  in  the  differential 
case,  and  it  has  an  extension  in  the  quasi-discrete  case  that  is  quite 
illuminating. 

We  denote  Fourier  Transforms  by  a  wavy  overhead  bar,  e.g.: 


oo 

c 

u(k,t)  =      u(x, 


t)e"ikx  dx  ,  (6) 


etc.   Then  (l)-(3)  reduce  directly  to: 


an     'V        '"b 

^  =  fv  -  ikg  h  ,  (7) 


|  =  -f  1  ,  (8) 


^  •   -tk  H  S  ,  (9) 


together  with  initial  conditions: 


o 


u(x, 


u  =  u(k,0)  =       u(x,0)e  lkx  dx  ,  (10) 


etc.   Since  (6)- (9)  is  a  coupled  set  of  constant  coefficient  ordinary 
differential  equations,  it  can  be  solved  by  the  usual  process  of  finding 
the  eigenvalues  and  eigenvectors  of  the  coefficient  matrix.   This  leads 
straightforwardly  to: 


'v    2  ivt    3  -ivt 
u  •  —  e    -  —  e 

v        v 


~   .,       if    ivt   if    -ivt 
v  =  ikg  a1  +  -~2   a_e *  a"\e 

v  v 


y        .  kH    ivt   kH    -ivt 

h  -  f  a  -  ~2  a2e    ^^ 

v  v 


> 


(11) 


where  v  is  given  by  (5),  and  the  a   are  picked  to  satisfy  the  initial 
conditions.   Observe  that  the  e     and  e      terms  both  represent  the 
transforms  of  transients  in  the  time  domain.   Solving  the  initial 
conditions  for  the  a   ,  collecting  like  terms  in  (11),  and  simplifying 


yields: 


f  v  ikg  h 

^/,    \    ^  o  o 

u(k,t)  =  u  cos  vt  H sin  vt 

o  v  v 


sin  vt  , 


-\  )k2gH   f2 

v(k,t)  =  —  u  sin  vt  +  < — r-  +  — r  cos  v 
v  o        J22 

V  V 


y..       ._     ikH  ^     ,        .        ikHf 

h(k,t)  = u  sin   vt — 

v        o  2 

v 


t(  vq  +  Q&  |l  -   cos   vt|hQ    ,  ) 
< 1  -   cos   vt>v     +  <  ■ 


f2   ^2gH 

_  +  B_  CQS  Vt)h 

2     2        I  o 

V       V 


Observe  that  (12)  immediately  yields  by  inspection  the  transform  of 
the  steady  state  ("balanced")  solution: 


us 

(k) 

- 

0 

vs 

(k) 

= 

k2gH 
2 

V 

'V 
V 

o 

+ 

ikef 

2 

V 

hs 

(k). 

:    m 

ikHf 
2 

V 

o 

+ 

2  h, 

2 

£    ^    f  (ike  £  ^  ) 

h  =  v  +  — -  <— f1  h  -  v  } 

o    o    2  (  f   o  oj 


>(13) 


^    H    f2 
h  +  -7   *  — ~  ik  » 
of    2      f   o 


{fK-W 


These  can  be  immediately  inverted  by  the  convolution  theorem,  and  the  obser- 
vation that 


00 

J 


-|x|/X  -ikx 
e       e     dx 


2A      2X  2 

2  2     2 

1+kT    v 


to  yield: 


"> 


us(x)  =  0 


vs(x)  =  v(x,0)  + 


hs(x)  =  h(x,0) 


00 


-  x-s  /A 


f  9x 


(s,0)  -  v(s,0)  [  ds 


V(14) 


2 

2A  f 


CO 

J 


-  x-s  /A 


f  8x 


(s,0)  -  v(s,0)  >  ds 


J 


Solving  equations  (l)-(3)  by  this  approach  allows  us  to  easily 
observe  several  elementary  results: 

(1)  Since  u  (k,0)   does  not  contribute  to  any  of  the  steady  state 
solutions,  any  "imbalance"  in  u(x,0)   does  not  affect  the  final  balanced 
state. 

(2)  The  term  in  (14) 


£  ill 
f  3x 


(s,0)  +  v(s,0)>  , 


the  initial  value  of  the  quantity  usually  referred  to  as  the  ageostrophic 
wind,  can  be  considered  as  the  measure  of  "imbalance"  in  the  initial  state 
that  is  not  dissipated  at  the  final  state. 

(3)  At  steady  state,  the  effects  of  this  initial  imbalance  contribute 
most  strongly  in  the  immediate  neighborhood  of  the  imbalance  and  die  off 
exponentially  in  space  away  from  it. 

Observe  also  that  the  transient  part  of  (12)  can  be  written: 


u  (k,t)  p  u  cos 
T  O 


f  K    ike  y 

vt  +  —  <v r°-  h  >  sin  vt 

v  )  o    f   o 


^  f  %  f  \  ^ 

vT(k,t)  =  -  -  uQsin  vt  +  — 


-^h 
f   o 


cos  vt 


£  ..       .  ikH  ^   .    t   ikHf  )^         ikg  y 

h,-,(k,t)  = u  sin  vt r-  <v 7s-  h  >  cos  vt 

T  v   o  2  )  o       o 


Wl5) 


Explicit  inverses  to  these  transforms  are  expressible  in  terms  of 

(l  I  2  2   2  2  J 

convolutions  involving  J  l"~v^  ft  -  x  /  .   However,  a  simplified  view  of 

the  asymptotic  behavior  of  the  transients  is  possible  by  using  the  method 

of  stationary  phase.   Let 

d(x,t)  -  v(x,t)  -*||  (x,t)  . 


Then 


t^Cx.t)  • 


h  f  **• 


0)cos  vt  e   dx  + 


hj    f*<k.?)rt 


ikx 
in  vt  e   dk  .  (16) 


Thus,  using  the  method  of  stationary  phase  (Appendix  1),  we  can  show,  for 
fixed  x  as  t  -*■  °°  : 


"t 


(x,t)   ^ 


2   2 

Xf   t 


,    ,2f22     2.3/2 
2tt(A    f   t   -x   ) 


< 


A2f2t2-x2 


cos 


k  iT2f  2t2-x2  + 


,2.2   2     2 

X    f   t   -x 


Xft 


Jx2f2t2-> 


,  o 


f  1      ^2-2   2     2 
sin  |  y  V  X   f   t  -x     +4', 


>-.        (17) 


where     <(>,       and     tk      are  slowly  varying  phases,      when      |x|    <   <   Xft    ,    this 


is  more   conveniently  approximated  as 


UT 


(x,t)   ^ 


2tt   X^ft 


'-fe'0 


cos(f  t+(j> ,  ) 


2C 


X2ft 


,  o 


sin(ft+i|; ,)  Y    . 


(18) 


Thus  the  decay  to  steady  state  of  u  at  a  fixed  point  x  is  governed 

by  two  factors: 

-1/2 

(a)  A  t      decay,  which  can  be  interpreted  as  the  effect  due  to 

the  dispersive  nature  of  the  process,  and 

(b)  An  additional  possible  decay,  depending  on  the  degree  to  which 
the  initial  imbalances  were  distributed  in  the  long  wave  numbers,  since 


|u  |-~ —  ,  0  |   =  Aim 


t-x»  ■   \X  ft 


t-*» 


(  u(s,0)e  i( 


-i(x/Xzft)s 


ds 


en 

"J 


u(s,0)ds  ,  etc. 


Similar  analysis  of     v   (x, t)      and     hT(x,t)      yield 


vT(x,t)   «b 


-i  1 
2 


2tt   A   ft 


«- 


a 


X 


A2ft 


,    0 


uf-f-,    0 

X   ft 


cos(ft+ik  )    /-   , 


sin(ft+<J>   ) 


(19) 


and 


hT(x,t)   % 


-i  1 


2tt   A   ft 


«- 


Hx 


2   2 
Aft 


rl 

u 


x2ft 


•)l 


sin(ft+<j>   ) 

(ft 


Hx 


2   2 
AZfZt 


31 


2       ■ 
A   ft 


cos(f  t+Ip    ) 


(20) 


Note   the  one  difference  in     h_(x, t).      Due  to   the  presence  of  the      (Ik) 

term  in     h    (k,t),    the  decay  due   to  dispersion  of     h    (x,t)      proceeds   as 

-3/2  -1/2 

t  ,    rather  than     t 

In  summary  then,  we  have  shown  here  that  the  differential  formulation 

of  the  dispersive  wave  model  (l)-(3)  is  most  easily  solved  by  Fourier 

transforming  to  a  system  of  ordinary  differential  equations.   This  solution 

shows  the  model  always  tends  to  a  steady  state,  whose  difference  from  the 

initial  state  is  determined  solely  by  the  initial  ageostrophic  wind  field. 

The  contributions  from  regions  where  the  ageostrophic  field  is  initally 

non-zero  die  out  in  the  steady  state  exponentially  with  distance.   Lastly, 

both  the  initial  u  and  ageostrophic  wind  fields  contribute  to  transients 

-1/2  -3/2 

that  die  out  in  time  as   t  '     for  the  velocities,  and  t  '     for  the 

free  surface  height. 


10 


3.    Second  Order  Centered  Finite  Differences  -  The  Quasi-Discrete  Adjust- 
ment 

In  this  section  we  present  an  analytic  treatment  of  the  most  common 
difference  scheme  used  for  (l)-(3),  and  investigate  the  resulting  effect 
on  the  steady-state  and  transient  behavior  discussed  above. 

Consider  the  continuous,  quasi-discrete,  second  order  centered-leap- 
frog  formulation  of  (l)-(3) 


—  =fv(x,t)   2Ax 


f-[h(x4-Ax,t)  -  h(x-Ax,t)] 


av 

3t 


3h 


=  -f  u(x,t) 


H 


8t     2Ax 


[u(x+Ax,t)  -  u(x-Ax, t)] 


>.(21) 


If  we  Fourier  Transform  this  system,  noting 


00 


u(x+Ax,t)e  dx 


ikAx  ^,        . 
e  u(x,t), 


we  have: 


"\ 


dt 


dt 


dh 
dt 


f  v  — ~°-  sin  kAx  h 

Ax 


a. 
-    f   u 


iH  ^ 

- —  sin  kAx  u 

Ax 


y  .(22) 
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Observe  this  system  is  identical  to  (7)- (9)  with  k  replaced  by: 


sin  kAx 
Ax 


(23) 


Thus  the  solution  to  (22)  becomes  identical  to  (12),  except  k  is  re- 
placed by  a    ,  and  v  by 


v  =  1  + 


(fe) 


4       21  A   I1/2 

sin  kAx 


(24) 


Observe  that  the  sinusoidal  terms  in  (12)  continue  to  represent  tran- 
sients  when  v  is  replaced  by  v  and  k  by  a  .   Thus  the  quasi  discrete 
case  will  tend  toward  a  steady  state  whose  transform  is: 


Vk) 


vs00 


0 


v  + 
o 


_J }J££  ft 

1+A  a     ' 


v 
o    o 


£  /,  N   £    H   1 

hs(k)  =ho  +  f  7TT2 

1+X  a 


lag  %    'x, 
la-;— r6-  h  -  v 
f   o    o 


>  (25) 


It  is  easily  shown  that 


i££  ft  -  v  -#")*  h(x+AxtQ)-h(x-AxlO) 
f   o    o    if 


2Ax 


-  v(x,0)[  .     (26) 


We  shall  let 


i(x,t)  =  ,<x.t>  _|MxH-Ax,t)-h(x-Ax,t) 

1  ^.^i  A 


(27) 


Thus  if  we  can  invert 


ia 


2  2 

1+A  a 


and 


1+A  a 


12 


the  steady  state   solution  would  be  available  by   convolution.      These   are 
inverted  in  Appendix  2,   where  we  show 

00  —6  |  x  |  /X  °° 

1         f  1 ikx  j  Ax  e \  A        . ,      0    .    N 

7T~       I       T~o   e  dx  =      , — —        /  5(x-2nAx), 

271     J       l+X2o2  /,21/A    ,2  ^ 


(Ax)  n 


S3— CO 


and 


00                                                       —6  |  x  |  /X  °° 

1         I  ia         ikx    ,,        -Ax  e \  A      x/      /0     -s.    N 


where 


•-fc-^r)   •  (28) 


and     0  < J   <   1     with  equality  only   for     Ax  =   0. 
Thus  we   can  arrive  at 


/°°    (     —3  |  s  I  / X  co 

)  077772    ^ 


v    (x)   =  v(x,0)   -  Ax  >  6(s-2nAx)l>d(x-s,0)ds 


„<         -6|2nAx|/A^ 
=  v(x,0) \       <>e  d(x-2nAx,0)(2Ax)>,    (29) 

2VA2+(Ax)2     n=-~ 


and  similarly 


„  V-^     (  -e|(2n-l)Ax|/A/N  ) 

hs(x)   =  h(x,0)  + ,       y        )e  d(x-(2n-l)Ax,0)(2Ax)V   . 

2AVX2+(Ax)2f     n=—   ( 

(30) 
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Clearly  (29)- (30)  tend  toward  the  corresponding  integral  forms  (14)  as 
Ax  ■*■  0. 

These  expressions  can  now  be  examined  versus  (14) .   The  conclusions 
that  we  can  draw  are  that  conversion  of  (l)-(3)  to  centered,  second  order, 
quasi-discrete  form  results  in: 

(1)  The  measure  of  imbalance,  the  initial  ageostrophic  wind,  is 
converted  to  a  finite  differenced  measure. 

(2)  Although  at  steady-state  the  effects  of  an  initial  disturbance 
die  off  exponentially  away  from  its  original  neighborhood,  the  rate  of 
decay  is  less  than  in  the  differential  case.   Furthermore,  the  rate  of 
decay  decreases  as   (Ax/ A)   increases. 

(3)  For  the  steady  states,  only  the  values  of  the  ageostrophic  wind 
at  alternating  points  are  considered. 

One  interpretation  of  why  the  sums  in  (29)- (30)  involves  only  alter- 
nating points  may  be  that,  when  elimination  is  tried  on  (21),  one  ends 
with  the  equation  for  u  involving  only  alternating  points: 


2 

—  u(x,t)  =  -f  u(x,t)  + — *sL_  iu(x+2Ax,t)-2u(x,t)+u(x-2Ax,t)>  . 

4(Ax)2  »  t         I 


A  complete  stationary  phase  analysis  of  the  transient  solution  in  the 
quasi-discrete  case: 


uT(k,t)  =  uQcos 


f  * 


vt  +  *i 


ice  ^ 
r  h  >  sin  vt 

v  (  °    f 


*      f   ;  *Vi      io"  2  ^ 

v^dtjt)  =  —  u  sin  vt  +  -7TX  <  v r*  h  >  cos  vt 

T       m  o        .,2)0    f   o 


W32) 


% 


hT(k,t)  = 


iaH  ^   .   -    iaHf  ~    ikg 
—  u  sin  vt rrr"  <  v 7^ 


K? 


J 
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is  algebraically  extremely  complicated,  although  quite  straightforward. 
However,  the  salient  features  are  relatively  easily  treated,  and  yield  the 
most  significant  results  on  the  transient  behavior.   Therefore  we  only 
present  an  outline  of  the  details. 

In  computing  the  inversion  integrals  for  (32),  terms  of  the  form 


oo 

f   A(k)e1(kx±0t) 


dk 


arise.   The  stationary  phase  points  of  these  integrals  arise  as  the 
solutions  of 


.    ,2,_     .  .   sin  kAx 
f  A  t  cos  kAx 


4>'(k)  =  x  ±     ,      ■  ■       Ax —  «  0  .       (33) 


1  + 


Ux")2sin2kAx 


This  expression  can  be  simplified  by  adding  (-x)  to  both  sides,  squaring 

2  2  2 

and  writing   cos  kAx  as   (1-sin  kAx),   to  yield  a  quadratic  in  sin  kAx. 

The  quadratic  formula  then  yields  as  the  points  of  stationary  phase  the 

solutions  of 

s1Aax-   aVt2-xV/(AW-*y-4fVxW_  ,      (34) 

2\    £    t 


It  is  then  easily  shown  that  as   t-*»  ,  these  solutions  closely  approximate: 

2          •  x2 

sin  kAx  =1  (35a) 

Aft 


and 


2,        2  -2 

.       2.     .  •  X      (AX)       A  /or,    \ 

sin  kAx  =  — T22 — 2 —     '  (35b) 

(A    f   t  -x  ) 
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Clearly,  for  large   t  ,  (35b)  yields  stationary  points  near  both 

.      ±  x  A"1 
k  = 


/2f22  2 
v\   f  t  -x 

and 


k-  ±7-. 

Ax 


(All  other  solutions  of  (35b)  are  beyond  the  Nyquist  cut-off.)   The  first 
points  are  slight  variants  of  the  stationary  points  for  the  differential 
case.   The  points  near  ±(tt/Ax)   arise  solely  from  the  discretization,  not 
the  physics  of  the  problem.   However,  since  these  points  give  a  behavior 
of   ((sin  kAx)/Ax)   identical  to  that  of  the  stationary  points  near  k  =  0, 
it  is  easily  shown  that  they  contribute  computational  transients  with 
precisely  the  same  asymptotic  behavior  as  the  physical  transients. 

Note,  before  we  consider  the  effect  of  terms  introduced  by  the  first 
solution,  (35a) ,  that  a  necessary  and  sufficient  condition  for  the 
stationary  phase  points  to  be  on  the  real  axis  is  that  the  quantity  under 
the  radical  sign  in  (34)  be  positive.   After  some  manipulation,  this 
condition  reduces  to 

lXl       (  /    2~~ 2       ) 
>   f  <  /(Ax)  +X   -  Ax  \  ,   t  >  0. 

It  is  not  coincidental  that  the  quantity  on  the  right  hand  side  of  the 
inequality  is  precisely  the  group  propagation  velocity  for  this  quasi- 
discrete  case. 

Referring  again  to  the  contribution  from  the  stationary  phase  points 
satisfying  (35a),  observe  that  we  can  easily  show  from  (33)  that 
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•TOO 


2 
±  A  ft 


1+(t—  )  sin  kAx 


(kf 


•,-(!+(— )   sin  kAx) sin  kAx  +  cos  kAx<- 


and  so,  near  these  stationary  points 


4»"(k)  ^ 


±A  ft  Ax 


/(Ax)2+X2 


Furthermore,  note  that  near  these  same  points 


/ 


2. ,2 


x;    /(Ax^+A 


Ax 


and    a   'v 


±1 
Ax 


Thus  following  again  the  argument  of  (3.7.5)  in  Miles  [5],  we  see  that  the 
points  of  stationary  phase  which  arise  from  the  solution  to  the  first  lead 
to  transients  whose  amplitudes,  asymptotically,  go  as  follows: 


u  cos   vt  -> 
o 


/ 


(Ax)2+A2 


_2tt   A   ft  Ax_ 


1/2 


*(*2S.°) 


f    U 
—  <v  - 


log,   £ 


o        f       o 


sin   vt-> 


Ax 


ll/2 


_2TTA2ft/(Ax)Z+A2_ 


^°)| 


f  ^      .      % 
—  u  sxn   vt 

v      ° 


Ax 


-il/2 


2TTA2ft/(Ax)2+ 


»(*&•<>) 


L.  h  -  ^  ft    i  cos  vt- 

,,2        o        f        o 


(Ax) 


1  1/2|* 


iaH  ^      . 
— —  u   sin   vt  -> 
o 


V 


_2ttA    ft 


H2f 


((Ax)2+A2) 


3/2 


^IS'0 


nl/2. 


_2TTA2tAx/(Ax)2+A2_ 


»'*^.o 
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iaHf  j~   lag  £  ) 


2  To   f     'C0S  " 


H2f  4x       11/2|* 
I  a 


°         o  ^2>.  //a  n2,,2\3/2 
_2ttA  t  ((Ax)  +A  J   . 


(**.o)| 


(Note  the  evaluations  are  at  ±(tt/2Ax)   since,  as  t-*»  ,  these  are  the 
only  solutions  of  (35a)  that  also  satisfy  the  Nyquist  limit.) 

Viewing  the  above,  it  is  now  clear  that  the  additional  stationary 
points  which  arise  in  this  quasi-discrete  case,  and  which  tend  toward 

±(tt/2Ax)   as  t-*30  cause  two  noticeable  effects  on  the  transient  behavior: 

-1/2 

(1)  All  of  these  transients  now  die  off  as  t      due  to  dispersion. 

Comparing  this  to  the  results  in  the  differential  case,  we  see  that  these 
transients  are  at  least  as  persistent  as  the  differential  transients,  and 

for  h(x,t)  more  so,  since  the  differential  transients  in  h(x,t)   die 

-3/2 
out  as  t 

(2)  Two  of  these  transients  (the  first  and  fifth)  are  somewhat  ill 
behaved  as   (Ax)  ■*  0  .   Assuming  u(x,0)  has  only  finite  power,  then  the 
Ax  term  in  the  denominator  should  be  controlled  by  tail-off  of 

u(±  -rr—  ,  0)   as   (Ax)  ■*-  0  ,  however,  these  terms  are  virtually  certain  to 
be  the  slowest  decaying  for  small  Ax  . 

(We  note  that  (35b)  also  causes  stationary  points  to  arise,  that  tend 
toward  ±(tt/Ax)   as   t-*»  ,  however,  the  transients  from  these  points  do 
not  have  an  asymptotic  dispersion  decay  that  depends  on   (Ax) ,  and  decay 
at  the  same  rate  as  the  differential  transients.)  Although  we  shall  not 
show  it  analytically,  we  suspect  the  additional  stationary  points  arise  in 
the  quasi-discrete  case  from  two  causes,  the  "folding"  in  temporal 
frequency  that  occurs  at  k  =  ±(tt/2Ax)  ,  and  the  high  frequency  cut  off 
at  k  =  ±(tt/Ax). 
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4.   Conclusion 

In  this  paper  we  have  examined  the  effect  of  second  order  centered 
spatial  discretization  on  a  dispersive  wave  equation.   We  have  shown  that 
the  methods  of  Fourier  Transforms,  and  in  particular  the  method  of  station- 
ary phase,  are  quite  useful  in  such  investigations.   In  both  the  differen- 
tial and  quasi-discrete  cases  we  have  provided  closed  form  expressions  for 
the  steady-state  solutions.   These  expressions  show  that  the  contribution 
from  any  point  in  the  initial  state  to  the  final  state  decays  exponentially 
with  distance  from  that  point  in  both  cases,  however  the  rate  of  decay  is 
decreased  by  discretization.   The  transients  in  both  cases  have  been 
analyzed  by  the  method  of  stationary  phase.   This  analysis  shows  that  the 
discretization  introduces  stationary  phase  points  that  have  no  counterpart 
in  the  differential  case,  and,  furthermore,  these  points  contribute 
transients  that  decay  no  faster  than  the  differential  ones,  and  in  one 
instance,  the  discrete  transient  will  dominate  the  differential  transient. 
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Appendix  1 


Consider  the  asymptotic  behavior  of 


oo 

2u  J 


A(k)cos  vt  eikxdk  =  i|)(x,t)  ,  (1-1) 


/   22 
where  v  -   f/l+X  k  .   This  integral  can  be  decomposed  into 


:x,t)  -AJJ   AC^e^l^dk  +   r  A(k)ei( 


where 


<J) ..  (k)  =  kx  +  vt,    and    <j>~(k)  =  kx  -  vt  . 


We  can  now  determine  the  asymptotic  behavior  as   t-*»  for  fixed  x 
of  this  term  by  using  the  method  of  stationary  phase.   Let  k-   and  k~ 


be  defined  by 


But 


4'l(kl)  =  Cf,2(k2)  =  °  ' 


<|>1(k)  =  x  + 


and  so  we  find  as  the  point  of  stationary  phase  for  the  first  integral: 


k. xA      ,  (1-2) 

which  is  on  the  path  of  integration,  for  t  >  x/Af  .   Similarly, 
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(f)'(k9)  =  0  yields  as  the  stationary  point  for  the  second  integral, 

, x/A 

k0  -  -  ,  (1-3) 

1        riA   2  2 

/A  f  t  -x 

which  is  also  on  the  path  of  integration  for  t  >  x/Af. 
Observe  that 


1      /?    2    2 2 
^(kj   =  -<|>0(k0)   =  t  A   f   t  -xz      , 


lv~l'  T2V    2'        A 

and   that 


,2.  M2*2.2      2n3/2 

^(kx)  =  -^(k2)  =      x  ft     =  (A  f  ^X  } —    >  0 


A-- 


'1+A   k. 


(Note      [♦JCk^l    «    U*^(k2)  |    =   ♦J0c1)    .    for      Aft   >    |x|.) 
Thus  we  have,   by    (3.7.5)    on  p.    51  of  Miles    [5]» 


2"*4fe)  K)e  " 


1/0/  ,-r    !    A2^2      2    ,    TT 

1/2  4  i[  —  /A   f   t  -x     +  -j 


+  A(k2)    e  X 


.r    1    /2,2   2 2        fr 
-i[  T  /A   f   t  -x     +  T  1 


/  o    o     W2    /  •  r    1    Z2.2   2      2        tt    , 

1  /         2tt  Af  2t2     \         (,.    .      i[yAft-x    +-] 

=  2\[A2f2t2-x2]3/2y    r   2  e 

-i[i/2f2t2-x2+|]l 
+  A(k2)    e  )    . 

But  note  that  if  A(k)   is  the  transform  of  a  purely  real  valued 
function, 
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■ 


a(x)eikxdk  =  A*(k)    . 


Thus,  we  can  write: 


ijj  ^ 


2   2 
Xf  -XT 


-,1/2 


.    ,,2,2   2     2.3/2 
2tt(X   f  t  -x  ) 


|  Re  A(k2)cos[  y  A2f2t2-x2  +  -| 


+  Im  A(k2)sin[  y  A2f2t2-x2  +  f  ]|   » 


(1-4) 


where     k~     is  given  by   (1-3).      Observe  for     Xft  >   >    |x|    ,    the  lead  term 


acts   as 


_2ttX   ft  _ 


1/2 


Similarly  we  can  show 


oo 

a.  r 

2tt      J 


A(k)sin  vt   e        dk 


'Xj 


1£2   2 
Xf   t 


-.1/2 


2,(X2f2t2-x2)3/2 


{t,      a  /i    \    •    r    1    /2.2   2     2       tt    . 
Re  A(k  )sxn[  -  Af   t  -x     +  j  J 


-   Im 


a  n     \         r    1   7(2.2^2      2        tt    n  ( 
A(k2)cos[  y  /X   f   t  -x     +  -r   ]  >  . 


(1-5) 


Lastly,    since   the  trigonetric  functions   in  both  expressions  have   the 
same  frequency,    and  since 


/[Re  A(k)]2  +    [Im  A(k)]2  =    |a(1c)  | 
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we  see 


00 

h  J 


A(k)cos   vt   e        dx  'v 


2   2 
Af   t 


.1/2 


„    ,, 2-2,2     2,3/2 

2 it  (A   f   t  -x  ) 


|A(k9)|cos[Y  A2f2t2-x2  +  f,] 


"2/,v-~olA 


?k" 


^k  =  |  -   tan  1[Im  A(k2)/Re  A(k2> ] 


=  j  -   arg(A(k2)) 


and 


oo 


A(k)sin  vt  e       dx  ^ 


,,2   2 
Af    t 


-a/2 


,,2.2,2     2.3/2 
2ir(A   f   t  -x  ) 


|A(k9)|sin[^  A2f2t2-x2  +  •.  ] 
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Appendix  2 


Consider 

00 

ikx 


1     f        e    dx     .         sin  k  Ax         ,_  ,,. 

2^      77T2~2  "  Where  °   =    Ax  (2~l) 

J       1  +  X  a 


The  denominator, 


*(k)  =  1  +  (jAZ   sin2k  Ax  (2-2) 


is  an  entire  function  of  k  in  the  complex  plane.   The  zeros  of  ^(k)   are 
solutions  of: 

sin  kAx  =  +  i  —  , 

or,  letting  k  «  k  +  i  k. 


sin(k  Ax  +  i  k.Ax) 
r       i 


=  sin(k  Ax)cosh(k.Ax)+i  cos(k  Ax)sinh(k.Ax) 
r        i  r        1 


■  ■  ■(?)  • 


Thus 


sin(k  Ax)  cosh (k.Ax)  =  0 
r         i 


cos 


(krAx)  sinh(k  Ax)  =  ±(p)  . 


Thus 


k  .±g,   k.  -  ±  f  sinh'Vf )   .  (2-3) 

r     Ax     i     Ax      yA  J 


as  opposed  to  the  poles  at  k  =  0,  k.  =  ±  1/A  for  the  continuous  case. 
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Thus  all  zeros  of  the  denominator  lie  off  the  real  axis.   In  fact,  they  are 
specifically  distributed  as  shown: 


a 

9 

<s 

1)          » 

0 

-3TT 

-2* 

-It 

TT 

AT 

a 

» 

(fe 

(3)              <b 

® 

Observe  further  that  at  its  zeros, 


*'  (k)  =  2  - —  sin  kAx  cos  kAx 

Ax 


=  ±2 


Ax 


(i  f )  coshlW1^ | 


=  ±2i/x2+(Ax)2 


(2-4) 


where  the  positive  sign  holds  for  k.  >  0,  and  the  negative  one  for  k  <0. 
Thus  the  zeros  of  iKk)   are  simple  poles  of  the  integrand. 
Consider  (2-1)  for  x  >  0.   (Observe  that  (2-1)  is  an 
even  function  of  x.)   To  insure  Re(ikx)  <0  we  will  close  the  contour 
in  Im(k)  >  0  half  plane.   As  long  as  we  route  the  contour  to  avoid  all 
the  poles  there,  the  integrand  is  exponentially  decaying  as   |k|  ■*■  »  , 
Im(k)  >  0.   Thus,  we  invoke  Jordan's  lemma  to  yield,  for  x  >  0 


h      J  7~fi  =   1   2-J    Residues  I- 


ikx 


+  X2a2' 


..E 


ikx 


lm(k)>0  ^'(k) 


k  +  ik. 
r     i 
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2/x2+(Ax 


x=  E-H-+1-slnhT)) 

(Ax)  -°° 


gx  °°  .n-rrx 


2/x2+(Ax)2 


E- 


=T"     e  >        e  ,  (2-5) 


where 


If)  • 


i^sinl 


Observe  (2-5)  is  not  a  convergent  series  in  the  usual  sense.  However, 
viewed  as  a  generalized  function,  it  is  a  Fourier  Series  for  the  function, 
periodic  of  period  2Ax  ,  given  in  the  interval  -Ax  <  x  <  Ax   by 

°°      .nrrx 

S(x)  m       /    c  e       where  c  =1. 
Z— t       n  n 


Ax        .mtr 

~  Ax 
But  observe  c  =  ~zr—         I    S(r)e      dr  =  1,   all  n. 


:n=2^    j    S<r> 


-Ax 
Thus 

S(t)  =  2Ax6(t), 
and  so  its  periodic  extension  becomes 


£  6(X-: 


S(x)  =  2Ax    7   6(x-2nAx). 


Thus 
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-3  x 
ikx 


—   I    r— -z   dk  =   ■     e        /       6(x-2nAx) . 


V^ 


(Ax) 


Note,  if  we  view   3  as  a  function 


3(x)  =  —  sinh  x  , 


it  is  easily  shown: 

3(0)  =  1    and    3'  (x)  <  0  . 

Thus   0  <  3(x)  <_    1  with  equality  only  at   x  =  0. 
Since,  for  Imk  >  0  and  large 

ia        Ax         Ax 


, , ,  2  2     2  .        ,2  x  Imk 
1+A  a     A  sin  kx    a  e 


the  above  argument  can  be  essentially  repeated  for 


ia(k)e 


h    [  zrri dk  - i  E 


lm(k>0  ^'(k) 


-Bx  °°  .niTX 


2X\X2+(Ax)' 


E  <->n  * 


A        \    1    /    i\H  Ax  .      _ 

e  >       (-1)      e  ,      x  >    0 


-p  lx|        °°        .mr,    ,  .    N 

-l  a      V^    X^(3d-Ax) 

e 


2X\/X2+(Ax) 


E 


-3    x 
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(2-7) 


^  e     A  Y^     6(x-(2n-l)Ax)  (2-8) 

XVX2+(Ax)2 
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